
*
* Post MC analysis plot maker
* 5/20/24
* ----------------------------------------------------------------------------

cd "~/Dropbox/Benton-Jordan-Philips/final-JOP-replication"
* ----------------------------------------------------------------------------


* ----------------------------------------------------------------------------
* Experiment I
* ----------------------------------------------------------------------------
use "data/dgp_scenario1_BASEPRED.dta", clear 
* already all labeled so join in the next:
joinby sims using "data/dgp_scenario1_MEPRED.dta", unmatched(both)
* Rename;, they're ME not PBS:
rename pred_median_PBS pred_median_ME
rename pred_ll_PBS pred_ll_ME
rename pred_ul_PBS pred_ul_ME
foreach var of varlist archterm_mean garchterm_mean hetx_mean hetz_mean constant_mean archterm_sd garchterm_sd hetx_sd hetz_sd constant_sd pred_sd {
	rename `var' `var'_ME
}
drop _merge
joinby sims using "data/dgp_scenario1_RESIDPRED.dta", unmatched(both)
* These names got messed up, they're RESID not PBS:
rename pred_median_PBS pred_median_RESID
rename pred_ll_PBS pred_ll_RESID
rename pred_ul_PBS pred_ul_RESID
foreach var of varlist archterm_mean garchterm_mean hetx_mean hetz_mean constant_mean archterm_sd garchterm_sd hetx_sd hetz_sd constant_sd pred_sd {
	rename `var' `var'_RESID
}
drop _merge
joinby sims using "data/dgp_scenario1_PBSPRED.dta", unmatched(both)
foreach var of varlist archterm_mean garchterm_mean hetx_mean hetz_mean constant_mean archterm_sd garchterm_sd hetx_sd hetz_sd constant_sd pred_sd {
	rename `var' `var'_PBS
}
drop _merge
* now join in delta method (this only has pred, se, ul and ll)
joinby sims using "data/dgp_scenario1_DELTAPRED.dta", unmatched(both)
drop _merge


* -----------------------------------------------------------------------------
*		Table SM.1, Experiment I Rows
* -----------------------------------------------------------------------------
* for median prediction
foreach var of varlist pred_median_PBS pred_DELTA pred_median_RESID pred_median_ME  {
	gen rmse_`var' = (`var' - pred_median_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
} 

* SD:
foreach var of varlist pred_sd_PBS pred_sd_DELTA pred_sd_RESID pred_sd_ME  {
	gen rmse_`var' = (`var' - pred_sd_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
} 

* UL:
foreach var of varlist pred_ul_PBS pred_ul_DELTA pred_ul_RESID pred_ul_ME  {
	gen rmse_`var' = (`var' - pred_ul_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
}

* LL:
foreach var of varlist pred_ll_PBS pred_ll_DELTA pred_ll_RESID pred_ll_ME  {
	gen rmse_`var' = (`var' - pred_ll_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
}
* -----------------------------------------------------------------------------


* -----------------------------------------------------------------------------
*		Figure SM.7 (a)
* -----------------------------------------------------------------------------
su pred_median_BASE, meanonly
graph hbox pred_median_PBS pred_median_RESID pred_median_ME pred_DELTA, nooutsides ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Predicted Median Expected Values") ///
 note("Note: Vertical line is calculated actual median prediction. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen1_predictions-box-wdelta.pdf", as(pdf) replace
* -----------------------------------------------------------------------------

* -----------------------------------------------------------------------------
*		Figure SM.7 (b)
* -----------------------------------------------------------------------------
su pred_sd_BASE, meanonly
graph hbox pred_sd_PBS pred_sd_RESID pred_sd_ME pred_sd_DELTA, nooutsides ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Standard Deviations") ///
 note("Note: Vertical line is calculated actual standard deviation. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen1_sd-box-wdelta.pdf", as(pdf) replace
* -----------------------------------------------------------------------------

* -----------------------------------------------------------------------------
*		Figure SM.7 (c)
* -----------------------------------------------------------------------------
su pred_ul_BASE, meanonly
graph hbox pred_ul_PBS pred_ul_RESID pred_ul_ME pred_ul_DELTA, nooutsides ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Upper 95% CI") ///
 note("Note: Vertical line is calculated actual upper 95% CI. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen1_ul-box-wdelta.pdf", as(pdf) replace
* -----------------------------------------------------------------------------

* -----------------------------------------------------------------------------
*		Figure SM.7 (d)
* -----------------------------------------------------------------------------
su pred_ll_BASE, meanonly
graph hbox pred_ll_PBS pred_ll_RESID pred_ll_ME pred_ll_DELTA, nooutsides ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Lower 95% CI") ///
 note("Note: Vertical line is calculated actual lower 95% CI. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen1_ll-box-wdelta.pdf", as(pdf) replace
* -----------------------------------------------------------------------------






* -------------------------------------------------------------------------
* Experiment II
* -------------------------------------------------------------------------

* join it all together
use "data/dgp_scenario2_BASEPRED.dta", clear 
* already all labled so join in the next:
joinby sims using "data/dgp_scenario2_MEPRED.dta", unmatched(both)
* Rename;, they're ME not PBS:
rename pred_median_PBS pred_median_ME
rename pred_ll_PBS pred_ll_ME
rename pred_ul_PBS pred_ul_ME
foreach var of varlist archterm_mean garchterm_mean hetx_mean hetz_mean constant_mean archterm_sd garchterm_sd hetx_sd hetz_sd constant_sd pred_sd {
	rename `var' `var'_ME
}
drop _merge
joinby sims using "data/dgp_scenario2_RESIDPRED.dta", unmatched(both)
* Rename; they're RESID not PBS:
rename pred_median_PBS pred_median_RESID
rename pred_ll_PBS pred_ll_RESID
rename pred_ul_PBS pred_ul_RESID
foreach var of varlist archterm_mean garchterm_mean hetx_mean hetz_mean constant_mean archterm_sd garchterm_sd hetx_sd hetz_sd constant_sd pred_sd {
	rename `var' `var'_RESID
}
drop _merge
joinby sims using "data/dgp_scenario2_PBSPRED.dta", unmatched(both)
foreach var of varlist archterm_mean garchterm_mean hetx_mean hetz_mean constant_mean archterm_sd garchterm_sd hetx_sd hetz_sd constant_sd pred_sd {
	rename `var' `var'_PBS
}
drop _merge
* now join in delta method (this only has pred, se, ul and ll)
joinby sims using "data/dgp_scenario2_DELTAPRED.dta", unmatched(both)
drop _merge

* -----------------------------------------------------------------------------
* 	Table SM.1 (Experiment II rows)
* -----------------------------------------------------------------------------
* for median prediction
foreach var of varlist pred_median_PBS pred_DELTA pred_median_RESID pred_median_ME  {
	gen rmse_`var' = (`var' - pred_median_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
} 

* SD:
foreach var of varlist pred_sd_PBS pred_sd_DELTA pred_sd_RESID pred_sd_ME  {
	gen rmse_`var' = (`var' - pred_sd_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
} 

* UL:
foreach var of varlist pred_ul_PBS pred_ul_DELTA pred_ul_RESID pred_ul_ME  {
	gen rmse_`var' = (`var' - pred_ul_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
}

* LL:
foreach var of varlist pred_ll_PBS pred_ll_DELTA pred_ll_RESID pred_ll_ME  {
	gen rmse_`var' = (`var' - pred_ll_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
}

* -----------------------------------------------------------------------------
*		Figure SM.8 (a)
* -----------------------------------------------------------------------------
su pred_median_BASE, meanonly
graph hbox pred_median_PBS pred_median_RESID pred_median_ME pred_DELTA, nooutsides ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Predicted Median Expected Values") ///
 note("Note: Vertical line is calculated actual median prediction. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen2_predictions-box-wdelta.pdf", as(pdf) replace
* -----------------------------------------------------------------------------


* -----------------------------------------------------------------------------
*		Figure SM.8 (b)
* -----------------------------------------------------------------------------
su pred_sd_BASE, meanonly
graph hbox pred_sd_PBS pred_sd_RESID pred_sd_ME pred_sd_DELTA, nooutsides ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Standard Deviations") ///
 note("Note: Vertical line is calculated actual standard deviation. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen2_sd-box-wdelta.pdf", as(pdf) replace
* -----------------------------------------------------------------------------


* -----------------------------------------------------------------------------
*		Figure SM.8 (c)
* -----------------------------------------------------------------------------
su pred_ul_BASE, meanonly
graph hbox pred_ul_PBS pred_ul_RESID pred_ul_ME pred_ul_DELTA, nooutsides ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Upper 95% CI") ///
 note("Note: Vertical line is calculated actual upper 95% CI. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen2_ul-box-wdelta.pdf", as(pdf) replace
* -----------------------------------------------------------------------------


* -----------------------------------------------------------------------------
*		Figure SM.8 (d)
* -----------------------------------------------------------------------------
su pred_ll_BASE, meanonly
graph hbox pred_ll_PBS pred_ll_RESID pred_ll_ME pred_ll_DELTA, nooutsides ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Lower 95% CI") ///
 note("Note: Vertical line is calculated actual lower 95% CI. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen2_ll-box-wdelta.pdf", as(pdf) replace
* -----------------------------------------------------------------------------







* ----------------------------------------------------------------------------
* Experiment III
* ----------------------------------------------------------------------------
use "data/dgp_scenario1t_BASEPRED.dta", clear 
* already all lableed so join in the next:
joinby sims using "data/dgp_scenario1t_MEPRED.dta", unmatched(both)
* Rename; they're ME not PBS:
rename pred_median_PBS pred_median_ME
rename pred_ll_PBS pred_ll_ME
rename pred_ul_PBS pred_ul_ME
foreach var of varlist archterm_mean garchterm_mean hetx_mean hetz_mean constant_mean archterm_sd garchterm_sd hetx_sd hetz_sd constant_sd pred_sd {
	rename `var' `var'_ME
}
drop _merge
joinby sims using "data/dgp_scenario1t_RESIDPRED.dta", unmatched(both)
* Rename; they're RESID not PBS:
rename pred_median_PBS pred_median_RESID
rename pred_ll_PBS pred_ll_RESID
rename pred_ul_PBS pred_ul_RESID
foreach var of varlist archterm_mean garchterm_mean hetx_mean hetz_mean constant_mean archterm_sd garchterm_sd hetx_sd hetz_sd constant_sd pred_sd {
	rename `var' `var'_RESID
}
drop _merge
joinby sims using "data/dgp_scenario1t_PBSPRED.dta", unmatched(both)
foreach var of varlist archterm_mean garchterm_mean hetx_mean hetz_mean constant_mean archterm_sd garchterm_sd hetx_sd hetz_sd constant_sd pred_sd {
	rename `var' `var'_PBS
}
drop _merge
* now join in delta method (this only has pred, se, ul and ll)
joinby sims using "data/dgp_scenario1t_DELTAPRED.dta", unmatched(both)
drop _merge


* -----------------------------------------------------------------------------
* 		Table SM.1 (Experiment III rows)
* -----------------------------------------------------------------------------
* for median prediction
foreach var of varlist pred_median_PBS pred_DELTA pred_median_RESID pred_median_ME  {
	gen rmse_`var' = (`var' - pred_median_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
} 

* SD:
foreach var of varlist pred_sd_PBS pred_sd_DELTA pred_sd_RESID pred_sd_ME  {
	gen rmse_`var' = (`var' - pred_sd_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
} 

* UL:
foreach var of varlist pred_ul_PBS pred_ul_DELTA pred_ul_RESID pred_ul_ME  {
	gen rmse_`var' = (`var' - pred_ul_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
}

* LL:
foreach var of varlist pred_ll_PBS pred_ll_DELTA pred_ll_RESID pred_ll_ME  {
	gen rmse_`var' = (`var' - pred_ll_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
}


* ------------------------------------------------------------------------------
*		Figure SM.9 (a)
* ------------------------------------------------------------------------------
su pred_median_BASE, meanonly
graph hbox pred_median_PBS pred_median_RESID pred_median_ME pred_DELTA, nooutsides ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Predicted Median Expected Values") ///
 note("Note: Vertical line is calculated actual median prediction. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen1t_predictions-box-wdelta.pdf", as(pdf) replace
* ------------------------------------------------------------------------------

* ------------------------------------------------------------------------------
*		Figure SM.9 (b)
* ------------------------------------------------------------------------------
su pred_sd_BASE, meanonly
graph hbox pred_sd_PBS pred_sd_RESID pred_sd_ME pred_sd_DELTA, nooutsides yscale(range(400)) ylabel(0(100)400) ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Standard Deviations") ///
 note("Note: Vertical line is calculated actual standard deviation. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen1t_sd-box-wdelta.pdf", as(pdf) replace
* ------------------------------------------------------------------------------

* ------------------------------------------------------------------------------
*		Figure SM.9 (c)
* ------------------------------------------------------------------------------
su pred_ul_BASE, meanonly
graph hbox pred_ul_PBS pred_ul_RESID pred_ul_ME pred_ul_DELTA, nooutsides ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Upper 95% CI") ///
 note("Note: Vertical line is calculated actual upper 95% CI. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen1t_ul-box-wdelta.pdf", as(pdf) replace
* ------------------------------------------------------------------------------

* ------------------------------------------------------------------------------
*		Figure SM.9 (d)
* ------------------------------------------------------------------------------
su pred_ll_BASE, meanonly
graph hbox pred_ll_PBS pred_ll_RESID pred_ll_ME pred_ll_DELTA, nooutsides ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Lower 95% CI") ///
 note("Note: Vertical line is calculated actual lower 95% CI. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen1t_ll-box-wdelta.pdf", as(pdf) replace
* ------------------------------------------------------------------------------







* -------------------------------------------------------------------------
* Experiment IV
* -------------------------------------------------------------------------
use "data/dgp_scenario2t_BASEPRED.dta", clear 
* already all labled so join in the next:
joinby sims using "data/dgp_scenario2t_MEPRED.dta", unmatched(both)
* Rename; they're ME not PBS:
rename pred_median_PBS pred_median_ME
rename pred_ll_PBS pred_ll_ME
rename pred_ul_PBS pred_ul_ME
foreach var of varlist archterm_mean garchterm_mean hetx_mean hetz_mean constant_mean archterm_sd garchterm_sd hetx_sd hetz_sd constant_sd pred_sd {
	rename `var' `var'_ME
}
drop _merge
joinby sims using "data/dgp_scenario2t_RESIDPRED.dta", unmatched(both)
* Rename; they're RESID not PBS:
rename pred_median_PBS pred_median_RESID
rename pred_ll_PBS pred_ll_RESID
rename pred_ul_PBS pred_ul_RESID
foreach var of varlist archterm_mean garchterm_mean hetx_mean hetz_mean constant_mean archterm_sd garchterm_sd hetx_sd hetz_sd constant_sd pred_sd {
	rename `var' `var'_RESID
}
drop _merge
joinby sims using "data/dgp_scenario2t_PBSPRED.dta", unmatched(both)
foreach var of varlist archterm_mean garchterm_mean hetx_mean hetz_mean constant_mean archterm_sd garchterm_sd hetx_sd hetz_sd constant_sd pred_sd {
	rename `var' `var'_PBS
}
drop _merge
* now join in delta method (this only has pred, se, ul and ll)
joinby sims using "data/dgp_scenario2t_DELTAPRED.dta", unmatched(both)
drop _merge


* ------------------------------------------------------------------------------
*		Table SM.1 (Experiment IV rows)
* ------------------------------------------------------------------------------
* for median prediction
foreach var of varlist pred_median_PBS pred_DELTA pred_median_RESID pred_median_ME  {
	gen rmse_`var' = (`var' - pred_median_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
} 

* SD:
foreach var of varlist pred_sd_PBS pred_sd_DELTA pred_sd_RESID pred_sd_ME  {
	gen rmse_`var' = (`var' - pred_sd_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
} 

* UL:
foreach var of varlist pred_ul_PBS pred_ul_DELTA pred_ul_RESID pred_ul_ME  {
	gen rmse_`var' = (`var' - pred_ul_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
}

* LL:
foreach var of varlist pred_ll_PBS pred_ll_DELTA pred_ll_RESID pred_ll_ME  {
	gen rmse_`var' = (`var' - pred_ll_BASE)^2 
	qui su rmse_`var', d
	loc rmse = round(sqrt(r(mean)), 0.001)
	loc rmedse = round(sqrt(r(p50)), 0.001)
	di "--"
	di "Var is `var' with RMSE: `rmse' and RMedSE: `rmedse'"
}


* ------------------------------------------------------------------------------
*		Figure SM.10 (a)
* ------------------------------------------------------------------------------
su pred_median_BASE, meanonly
graph hbox pred_median_PBS pred_median_RESID pred_median_ME pred_DELTA, nooutsides ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Predicted Median Expected Values") ///
 note("Note: Vertical line is calculated actual median prediction. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen2t_predictions-box-wdelta.pdf", as(pdf) replace
* ------------------------------------------------------------------------------


* ------------------------------------------------------------------------------
*		Figure SM.10 (b)
* ------------------------------------------------------------------------------
su pred_sd_BASE, meanonly
graph hbox pred_sd_PBS pred_sd_RESID pred_sd_ME pred_sd_DELTA, nooutsides ///
 yscale(range(100)) ylabel(0(20)100) ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Standard Deviations") ///
 note("Note: Vertical line is calculated actual standard deviation. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen2t_sd-box-wdelta.pdf", as(pdf) replace
* ------------------------------------------------------------------------------


* ------------------------------------------------------------------------------
*		Figure SM.10 (c)
* ------------------------------------------------------------------------------
su pred_ul_BASE, meanonly
graph hbox pred_ul_PBS pred_ul_RESID pred_ul_ME pred_ul_DELTA, nooutsides ///
	yscale(range(200)) ylabel(0(50)200) ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Upper 95% CI") ///
 note("Note: Vertical line is calculated actual upper 95% CI. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen2t_ul-box-wdelta.pdf", as(pdf) replace
* ------------------------------------------------------------------------------


* ------------------------------------------------------------------------------
*		Figure SM.10 (d)
* ------------------------------------------------------------------------------
su pred_ll_BASE, meanonly
graph hbox pred_ll_PBS pred_ll_RESID pred_ll_ME pred_ll_DELTA, nooutsides ///
 yline(`r(mean)', lcolor(black) lpattern(solid) lwidth(medthick)) ///
 box(1, color(gray%75)) box(2, color(edkblue%75)) box(3, color(eltgreen%75)) ///
 box(4, color(cranberry%75)) ///
 title("Dist. of 500 Bootstrap Lower 95% CI") ///
 note("Note: Vertical line is calculated actual lower 95% CI. Extreme values omitted for clarity.") /// 
 legend(order(1 "Parametric" 2 "Residual" 3 "Maximum Entropy" 4 "Delta") rows(1)) ///
 ytitle("Prediction")
*graph export "scen2t_ll-box-wdelta.pdf", as(pdf) replace
* ------------------------------------------------------------------------------



